function x = SorSolver(lhs, rhs, numIter)
    w = 0.7;
    
    x = zeros(size(lhs, 2), 1);
    L = sparse(tril(lhs, -1));
    U = sparse(triu(lhs, 1));
    D = spdiags(spdiags(lhs, 0), 0, size(lhs, 1), size(lhs, 2));
    A1 = D + w * L;
    A2 = (w * U + (w - 1) * D);
    scaledRhs = w * rhs;
      
    for i = 1 : numIter        
        x = A1 \ (scaledRhs -  A2 * x);   
    end
end